home *** CD-ROM | disk | FTP | other *** search
/ Die Ultimative Software-P…i Collection 1996 & 1997 / Die Ultimative Software-Pakete CD-ROM fur Atari Collection 1996 & 1997.iso / g / gnu_c / gpplib22.zoo / libsrc / xrationa.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1993-11-30  |  7.3 KB  |  415 lines
  1. /* 
  2. Copyright (C) 1988 Free Software Foundation
  3.     written by Doug Lea (dl@rocky.oswego.edu)
  4.  
  5. This file is part of the GNU C++ Library.  This library is free
  6. software; you can redistribute it and/or modify it under the terms of
  7. the GNU Library General Public License as published by the Free
  8. Software Foundation; either version 2 of the License, or (at your
  9. option) any later version.  This library is distributed in the hope
  10. that it will be useful, but WITHOUT ANY WARRANTY; without even the
  11. implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  12. PURPOSE.  See the GNU Library General Public License for more details.
  13. You should have received a copy of the GNU Library General Public
  14. License along with this library; if not, write to the Free Software
  15. Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  16. */
  17.  
  18. #ifdef __GNUG__
  19. #pragma implementation
  20. #endif
  21. #include <xrationa.h>
  22. #include <std.h>
  23. #include <math.h>
  24. #include <builtin.h>
  25. #include <float.h>
  26.  
  27. void Rational::error(const char* msg) const
  28. {
  29.   (*lib_error_handler)("Rational", msg);
  30. }
  31.  
  32. static const Integer _Int_One(1);
  33.  
  34. void Rational::normalize()
  35. {
  36.   int s = sign(den);
  37.   if (s == 0)
  38.     error("Zero denominator.");
  39.   else if (s < 0)
  40.   {
  41.     den.negate();
  42.     num.negate();
  43.   }
  44.  
  45.   Integer g = gcd(num, den);
  46.   if (ucompare(g, _Int_One) != 0)
  47.   {
  48.     num /= g;
  49.     den /= g;
  50.   }
  51. }
  52.  
  53. void      add(const Rational& x, const Rational& y, Rational& r)
  54. {
  55.   if (&r != &x && &r != &y)
  56.   {
  57.     mul(x.num, y.den, r.num);
  58.     mul(x.den, y.num, r.den);
  59.     add(r.num, r.den, r.num);
  60.     mul(x.den, y.den, r.den);
  61.   }
  62.   else
  63.   {
  64.     Integer tmp;
  65.     mul(x.den, y.num, tmp);
  66.     mul(x.num, y.den, r.num);
  67.     add(r.num, tmp, r.num);
  68.     mul(x.den, y.den, r.den);
  69.   }
  70.   r.normalize();
  71. }
  72.  
  73. void      sub(const Rational& x, const Rational& y, Rational& r)
  74. {
  75.   if (&r != &x && &r != &y)
  76.   {
  77.     mul(x.num, y.den, r.num);
  78.     mul(x.den, y.num, r.den);
  79.     sub(r.num, r.den, r.num);
  80.     mul(x.den, y.den, r.den);
  81.   }
  82.   else
  83.   {
  84.     Integer tmp;
  85.     mul(x.den, y.num, tmp);
  86.     mul(x.num, y.den, r.num);
  87.     sub(r.num, tmp, r.num);
  88.     mul(x.den, y.den, r.den);
  89.   }
  90.   r.normalize();
  91. }
  92.  
  93. void      mul(const Rational& x, const Rational& y, Rational& r)
  94. {
  95.   mul(x.num, y.num, r.num);
  96.   mul(x.den, y.den, r.den);
  97.   r.normalize();
  98. }
  99.  
  100. void      div(const Rational& x, const Rational& y, Rational& r)
  101. {
  102.   if (&r != &x && &r != &y)
  103.   {
  104.     mul(x.num, y.den, r.num);
  105.     mul(x.den, y.num, r.den);
  106.   }
  107.   else
  108.   {
  109.     Integer tmp;
  110.     mul(x.num, y.den, tmp);
  111.     mul(y.num, x.den, r.den);
  112.     r.num = tmp;
  113.   }
  114.   r.normalize();
  115. }
  116.  
  117.  
  118.  
  119.  
  120. void Rational::invert()
  121. {
  122.   Integer tmp = num;  
  123.   num = den;  
  124.   den = tmp;  
  125.   int s = sign(den);
  126.   if (s == 0)
  127.     error("Zero denominator.");
  128.   else if (s < 0)
  129.   {
  130.     den.negate();
  131.     num.negate();
  132.   }
  133. }
  134.  
  135. int compare(const Rational& x, const Rational& y)
  136. {
  137.   int xsgn = sign(x.num);
  138.   int ysgn = sign(y.num);
  139.   int d = xsgn - ysgn;
  140.   if (d == 0 && xsgn != 0) d = compare(x.num * y.den, x.den * y.num);
  141.   return d;
  142. }
  143.  
  144. Rational::Rational(double x)
  145. {
  146.   num = 0;
  147.   den = 1;
  148.   if (x != 0.0)
  149.   {
  150.     int neg = x < 0;
  151.     if (neg)
  152.       x = -x;
  153.  
  154.     const long shift = 15;         // a safe shift per step
  155.     const double width = 32768.0;  // = 2^shift
  156.     const int maxiter = 20;        // ought not be necessary, but just in case,
  157.                                    // max 300 bits of precision
  158.     int expt;
  159.     double mantissa = frexp(x, &expt);
  160.     long exponent = expt;
  161.     double intpart;
  162.     int k = 0;
  163.     while (mantissa != 0.0 && k++ < maxiter)
  164.     {
  165.       mantissa *= width;
  166.       mantissa = modf(mantissa, &intpart);
  167.       num <<= shift;
  168.       num += (long)intpart;
  169.       exponent -= shift;
  170.     }
  171.     if (exponent > 0)
  172.       num <<= exponent;
  173.     else if (exponent < 0)
  174.       den <<= -exponent;
  175.     if (neg)
  176.       num.negate();
  177.   }
  178.   normalize();
  179. }
  180.  
  181.  
  182. Integer trunc(const Rational& x)
  183. {
  184.   return x.num / x.den ;
  185. }
  186.  
  187.  
  188. Rational pow(const Rational& x, const Integer& y)
  189. {
  190.   long yy = y.as_long();
  191.   return pow(x, yy);
  192. }               
  193.  
  194. #if defined(__GNUG__) && !defined(NO_NRV)
  195.  
  196. Rational operator - (const Rational& x) return r(x)
  197. {
  198.   r.negate();
  199. }
  200.  
  201. Rational abs(const Rational& x) return r(x)
  202. {
  203.   if (sign(r.num) < 0) r.negate();
  204. }
  205.  
  206.  
  207. Rational sqr(const Rational& x) return r
  208. {
  209.   mul(x.num, x.num, r.num);
  210.   mul(x.den, x.den, r.den);
  211.   r.normalize();
  212. }
  213.  
  214. Integer floor(const Rational& x) return q
  215. {
  216.   Integer r;
  217.   divide(x.num, x.den, q, r);
  218.   if (sign(x.num) < 0 && sign(r) != 0) --q;
  219. }
  220.  
  221. Integer ceil(const Rational& x) return q
  222. {
  223.   Integer  r;
  224.   divide(x.num, x.den, q, r);
  225.   if (sign(x.num) >= 0 && sign(r) != 0) ++q;
  226. }
  227.  
  228. Integer round(const Rational& x) return q
  229. {
  230.   Integer r;
  231.   divide(x.num, x.den, q, r);
  232.   r <<= 1;
  233.   if (ucompare(r, x.den) >= 0)
  234.   {
  235.     if (sign(x.num) >= 0)
  236.       ++q;
  237.     else
  238.       --q;
  239.   }
  240. }
  241.  
  242. // power: no need to normalize since num & den already relatively prime
  243.  
  244. Rational pow(const Rational& x, long y) return r
  245. {
  246.   if (y >= 0)
  247.   {
  248.     pow(x.num, y, r.num);
  249.     pow(x.den, y, r.den);
  250.   }
  251.   else
  252.   {
  253.     y = -y;
  254.     pow(x.num, y, r.den);
  255.     pow(x.den, y, r.num);
  256.     if (sign(r.den) < 0)
  257.     {
  258.       r.num.negate();
  259.       r.den.negate();
  260.     }
  261.   }
  262. }
  263.  
  264. #else
  265.  
  266. Rational operator - (const Rational& x) 
  267. {
  268.   Rational r(x); r.negate(); return r;
  269. }
  270.  
  271. Rational abs(const Rational& x) 
  272. {
  273.   Rational r(x);
  274.   if (sign(r.num) < 0) r.negate();
  275.   return r;
  276. }
  277.  
  278.  
  279. Rational sqr(const Rational& x)
  280. {
  281.   Rational r;
  282.   mul(x.num, x.num, r.num);
  283.   mul(x.den, x.den, r.den);
  284.   r.normalize();
  285.   return r;
  286. }
  287.  
  288. Integer floor(const Rational& x)
  289. {
  290.   Integer q;
  291.   Integer r;
  292.   divide(x.num, x.den, q, r);
  293.   if (sign(x.num) < 0 && sign(r) != 0) --q;
  294.   return q;
  295. }
  296.  
  297. Integer ceil(const Rational& x)
  298. {
  299.   Integer q;
  300.   Integer  r;
  301.   divide(x.num, x.den, q, r);
  302.   if (sign(x.num) >= 0 && sign(r) != 0) ++q;
  303.   return q;
  304. }
  305.  
  306. Integer round(const Rational& x) 
  307. {
  308.   Integer q;
  309.   Integer r;
  310.   divide(x.num, x.den, q, r);
  311.   r <<= 1;
  312.   if (ucompare(r, x.den) >= 0)
  313.   {
  314.     if (sign(x.num) >= 0)
  315.       ++q;
  316.     else
  317.       --q;
  318.   }
  319.   return q;
  320. }
  321.  
  322. Rational pow(const Rational& x, long y)
  323. {
  324.   Rational r;
  325.   if (y >= 0)
  326.   {
  327.     pow(x.num, y, r.num);
  328.     pow(x.den, y, r.den);
  329.   }
  330.   else
  331.   {
  332.     y = -y;
  333.     pow(x.num, y, r.den);
  334.     pow(x.den, y, r.num);
  335.     if (sign(r.den) < 0)
  336.     {
  337.       r.num.negate();
  338.       r.den.negate();
  339.     }
  340.   }
  341.   return r;
  342. }
  343.  
  344. #endif
  345.  
  346. ostream& operator << (ostream& s, const Rational& y)
  347. {
  348.   if (y.denominator() == 1L)
  349.     s << y.numerator();
  350.   else
  351.   {
  352.     s << y.numerator();
  353.     s << "/";
  354.     s << y.denominator();
  355.   }
  356.   return s;
  357. }
  358.  
  359. istream& operator >> (istream& s, Rational& y)
  360. {
  361. #ifdef _OLD_STREAMS
  362.   if (!s.good())
  363.   {
  364.     return s;
  365.   }
  366. #else
  367.   if (!s.ipfx(0))
  368.   {
  369.     s.clear(ios::failbit|s.rdstate()); // Redundant if using GNU iostreams.
  370.     return s;
  371.   }
  372. #endif
  373.   Integer n = 0;
  374.   Integer d = 1;
  375.   if (s >> n)
  376.   {
  377.     char ch = 0;
  378.     s.get(ch);
  379.     if (ch == '/')
  380.     {
  381.       s >> d;
  382.     }
  383.     else
  384.     {
  385.       s.putback(ch);
  386.     }
  387.   }
  388.   y = Rational(n, d);
  389.   return s;
  390. }
  391.  
  392. int Rational::OK() const
  393. {
  394.   int v = num.OK() && den.OK(); // have valid num and denom
  395.   if (v)
  396.     {
  397.       v &= sign(den) > 0;           // denominator positive;
  398.       v &=  ucompare(gcd(num, den), _Int_One) == 0; // relatively prime
  399.     }
  400.   if (!v) error("invariant failure");
  401.   return v;
  402. }
  403.  
  404. int
  405. Rational::fits_in_float() const
  406. {
  407.     return Rational (FLT_MIN) <= *this && *this <= Rational (FLT_MAX);
  408. }
  409.  
  410. int
  411. Rational::fits_in_double() const
  412. {
  413.     return Rational (DBL_MIN) <= *this && *this <= Rational (DBL_MAX);
  414. }
  415.